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& ! Abstract 
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A new algorithm for determining the eigenstates of n-dimensional billiards 

O ■ is presented. It is based on the application of the Cauchy theorem for the 

determination of the null space of the boundary overlap matrix. The method is 

free from the limitations associated with the shape of the billiard and could be 

applied even for nonconvex geometries where other algorithms face difficulties. 

I/" - ) ' Moreover it does not suffer from the existence of eigenvalue degeneracies which 

is another serious shortcoming of many methods. In the paper we apply the 

algorithm to a few simple cases where the analytical solutions exist. Numerical 

ON , solutions have been investigated for the case of annular billiard. 
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I. INTRODUCTION 

In spite of the fact that determining the spectrum of a cavity is certainly more than a 
century old problem, numerical algorithms do not seem to abound, and moreover, they do 
not seem always to solve the problem entirely. The references to known methods familiar to 
us, do not seem to be always the original ones. However, these references apparently cover 
all the main ideas used so far, with the exception of a few rather specialized algorithms. 
We shall mention here only those methods which have the potential to be used for a rather 
large class of billiard shapes and hopefully in more than two dimensions. 

The first type of algorithms amounts to determining the roots of a very high order 
determinant 0-^[. This is not a very pleasant numerical task, as such a determinant easily 
can vary by many orders of magnitude and the numerically null space is certainly another 
major problem, which can and apparently does lead to spurious states. The variant suggested 
in Ref. is somewhat worse in some respect, as it requires determination of complex roots 
of a complex valued determinant. Even though the spectrum of a cavity is real, when 
implemented in an obvious and natural way the method used in Ref. || leads to spurious 
states, the nature of which is mysterious. Moreover it has been shown recently that these 
methods have also problems for nonconvex billiard geometries ||. 

In Ref. J7||| another method was used, which amounts to determining relatively sharp 
minima of the norm of the wave function over the domain boundary. Such minima do not 
seem to be always well defined as there are some spurious minima as well, not so deep 
however. Moreover it is not entirely clear how to disentangle single and multiple degenerate 
states. The determination of all eigenstates requires multiple runs with randomly chosen 
conditions ||. The only way to decide whether all the eigenvalues have been found is by 
comparing the determined average level density with the Weyl formula. Even though this 
method seems to reveal all the eigenvalues in the end, it does appear to us to be a rather 
time consuming procedure, needing a "constant outside intervention and evaluation" of the 
results and thus requiring too much time. 

A significantly superior method, which apparently is free of many of the above drawbacks 
was recently suggested in Ref. [|10| . One of the major deficiencies of this method however is 



its limitation to the so called star-shaped domains. A star-shaped domain is one for which ft- 
r > 0, where n is the outward normal to the domain boundary and r is the position vector of 
a point on the domain boundary with respect to an arbitrary origin within the domain. The 
restriction to star-shaped domains excludes for example annular billiards, peanut-shaped 
billiards and so forth, which are extremely interesting for various applications. Besides this, 
a few rather small technical details, concerning its actual numerical implementation have 
been only slightly covered in Ref. [TO and thus the numerical limitations of this method 
are not entirely clear. But in any case, this method apparently is the most powerful one 
available in literature so far. 



Another method known to us is the constraint operator method presented in Ref. |iT 
Although it allows to determine all eigenvalues in the truncated space, it does not seem to 
be very accurate for some states, especially the high lying states. The reason could be that 
the algorithm requires several numerical integrations and diagonalizations. It is also not 
obvious whether it can actually be used for domains with holes for example. 

In this paper we are suggesting an entirely different method, which when correctly imple- 



mented should lead to the determination of the whole eigenspectrum in a given interval. The 
method might not necessarily be the fastest, as the method suggested in Ref. [I(J is perhaps 



faster. However, the method presented here should work for arbitrary shapes, unlike the 
one in Ref. [Kj. Its numerical implementation does not seem to suffer from any particular 



and/or severe drawbacks. 

II. ALGORITHM 

The solution of the Helmholtz equation 

A0(r) = -kl^r) (2.1) 

in the n-dimensional connected domain T> with the Dirichlet boundary conditions on the 
boundary B 

<f>{r)\ B = (2.2) 

can be represented as follows 

N 

^) = Ecnifer) (2.3) 

n=\ 

where functions tp n (k ,r) are solutions of the Helmholtz equation fl2.1|) , which nevertheless 



do not satisfy the boundary condition ( |2.2j ). Usually it is the most convenient to choose 
functions ipn in the form of plane waves: exp{ikok n ■ r) with condition: k n ■ k n = 1 (see 
however Section IIIC). In such a case the sum in ( |2.3| ) should be understood as an integral 
over all possible orientations of the unit vectors k n . Since in any numerical realization 
an integral is always discretized in one way or another, we shall use a discrete summation 
throughout most of this paper. (The interested reader should find no difficulty in turning this 
"approximate" representation into an "exact" one, if such a need arises.) The eigenvectors 
and eigenvalues will be determined from the condition that the boundary norm vanishes, 
namely that 



1 / N N 

- i dSw(r) Y. c* n 4>* n (ko,r)^ m (k ,r)c m = Y. c* n O nm (k )c m = 0, (2.4) 

^ JB n,m=l n,m=l 

where we have chosen arbitrarily a unit normalization of the diagonal matrix elements, as 
S = § B dSw(r). An arbitrary positively defined weight function w(r) > can be introduced 
as well and sometimes this flexibility of the formalism can be profitably used. In this paper we 
shall limit ourselves to a unit weight function w(r) = 1, unless otherwise noted. Obviously, 
for any finite domain T> the quantization condition Eq. ( j2.4|) is satisfied only for discrete 
values of fco- The determination of these eigenvalues and of their degeneracies is numerically 
the most difficult problem. For arbitrary values of k let us introduce the eigenvalues and 
eigenvectors of the boundary overlap matrix O(k) (BOM) 

0(k)C a = \ a (k)C a , (2.5) 



where C a is a column vector and the appropriate matrix to vector multiplication is implied. 
From the non-negativity of the boundary norm it follows that for real k these eigenvalues 
are non-negative 

\ a {k) > 0. (2.6) 

Only for an eigenvalue k$ of the Helmholtz equation with the Dirichlet boundary conditions 
one or more BOM eigenvalues vanish, \ a (k ) = 0. The number of such eigenvalues \ a (k ), 
which vanish simultaneously, is equal to the degeneracy of the corresponding eigenvalue of 
the original Helmholtz equation. Even though we did not rigorously proved it yet, it appears 
that in the neighborhood of an eigenvalue of the Helmholtz equation, those BOM eigenvalues 
which vanish, always do so quadratically only, i.e. 

X a (k) oc(k- k ) 2 + . . . (2.7) 

The basic idea behind our approach is to analytically continue the BOM O(k) into the 
complex fc-plane, namely to make k complex in the definition of the matrix elements 

O nm {k) = i / dS^Jk, r)ip m (k } r) (2.8) 

and to compute around a contour C in the complex /c-plane the integral 

m= f»£m. (2 . 9) 

Here X' a (k) is the derivative of the eigenvalue X a (k) with respect to k. This derivative can 
be easily calculated once the eigenvectors of the BOM are known as follows 

\' a (k) = C a (kyO'(k)C a (k), (2.10) 

where 0'(k) is the derivative of 0(k) with respect to k, C a (k) are normalized as usual 
C^(k)C a (k) = 1 and in all these relations obvious vector and matrix multiplications are 
implied. As an attentive reader would have remarked, we have divided the value of the 
integral by 4iii, instead of the naively expected 2iri, in order to take into account the double 
degeneracy of a root of X a {k), we have discussed above, see Rels. ( |2.6|) and Q2.71) . Also, when 
k is complex, when computing the BOM elements the 6ra-vectors are actually functions of 
k* , while the A;et-vectors are functions of k. Only by defining the BOM elements in this way 
O nm (k) is an analytical function, therefore a function of k only and not a function of k and 
k* as well. We shall consider here contours C in the semiplane Re (A;) > only. The integral 
in Rel. ( |2.9p can be performed analytically, as 

r Hk N \' (k\ 1 N 

We do not use this obvious result, as in any computer implementation the determination of 
the actual Riemannian sheet and of the change of the logarithm around the contour seems 



to be ambiguous. The numerical evaluation of the integral however, appears to be always 
straightforward to implement. 

M{C) is thus exactly equal to the number of eigenvalues (counting the number of degen- 
eracies as well), of the Helmholtz equation Eq. ([2.1|) with the Dirichlet boundary conditions 
Eq. ( p.2| ), on the segment of the real fc-axis enclosed by the contour C. Consequently, no 
eigenvalue of the original equation can thus be missed. This statement is strictly speaking 
correct only in the limit N — ■> oo. However, one can easily convince oneself that if N is 
suitable large, this holds true anyway. Note that M{C) takes only integer values and thus 
when N is larger than a certain value the N — ► oo limit is exactly attained (unless the unit 
vectors k n are distributed in a very peculiar manner of course) (see table [[]). 

One can introduce also other useful quantities. For example 

^=/ c 3££xi*" = £* (2 ' 12) 
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where the sum over /3 G C is over all eigenvalues enclosed by the contour C and n is in 
principle an arbitrary number, not necessarily an integer and positive (N.B. the origin is 
not encircled by the contour C) . If the contour C encloses only one eigenvalue then its value 
is equal to kp = Si(C). By determining S n {C) for n — 0, . . . , u, where v = So(C) = Af(C), 
one can easily set up a polynomial, whose roots are the eigenvalues enclosed inside the 
contour C. The functions S n (C) have similar convergence properties as N(C) discussed 
above. Obviously, one can compute in this way arbitrary functions of the eigenspectrum as 
well. In particular, it is useful to calculate Fourier components of the quantity S n {C)\ 

4(C «-> = { £ t t ^'V- = g^M~, (, 13) 

Tfl 

where t m = 27r——,m = —L, ...,L and Ak is the length of the interval on the real k— axis 

Ak 
enclosed by the contour C. Thus through the inverse Fourier transform one can obtain e.g. 

level density p or energy distributions g inside the contour C: 

1 L 



Pc{k) = 2LTT £ MC,t m )e' M -, 

rn=—L 

9c{k) = 2FTl ^ ^(C,t m )e- M -. (2.14) 



m=—L 



Here we put pc(k) = gc(k) — for k outside the contour. Since the numerical costs of 
calculations of the Fourier transform S n and S n are of the same order we can gain at almost 
no expense more precise information about the distribution of the functions S n (C) inside 
the integration contour. 

The practical implementation of this algorithm is rather straightforward. There are 
essentially a few relatively simple aspects one has to keep in mind: 

• One should divide the real /c-axis in not too long intervals to be enclosed by a complex 
contour C. When computing the BOM eigenvalues X a (k) along such a contour one 



should use always the same number of basis wave functions ip n . When the basis is 
increased, even though one might not gain in overall numerical accuracy, there is a 
side effect. Irrelevant eigenvalues X a (k) give a large contribution to the integrand in 
Rel. (|2]9|) and thus changing the number of wave functions ip n along the contour leads 
to erroneous results. 

The total number of basis wave functions should be chosen somewhat larger than the 
number of quantum states the boundary can accommodate in this energy range. For 
example, for a 2-dimensional boundary one should have of the order of Lk/n plane 
waves, where L is the length of the outer perimeter. 

As the method gives the exact number of the eigenstates in a given energy interval, one 
can easily narrow the interval so as to determine the exact location of any eigenstate 
and its degeneracy as well, along with the corresponding eigenvector of the Helmholtz 
equation. It is unavoidable to have parts of the contour C close to the real fc-axis. It is 
profitable however to choose the imaginary part of the contour C not too far and not 
to close to the real fc-axis. A symmetric rectangular contour, with two sides parallel 
to the real fc-axis and two sides normal to it seems like a most reasonable and flexible 
choice. Parts of such a contour can be used repeatedly in order to narrow down the 
position of the actual eigenvalues. 

The matrix elements of BOM can differ between each other by many orders of magni- 
tude if the basis is large. It can affect the numerical accuracy of the method. Therefore 
it is recommended in such a case to rescale the BOM matrix: O(k) = J r O{k)J r , where 
detF ^ 0. 

The imaginary part of the contour should be chosen at a distance of the order of the 
mean level separation. At closer distances the integrand changes too rapidly with k, 
while at larger distances from the real A;-axis a significantly increased and unnecessary 
numerical accuracy might be required to compute the integrand. 

The fact that M{C) is integer valued makes its computation somewhat easy, as a 
relatively low numerical accuracy can be used however. 

The Fourier transform calculated in the finite interval contains spurious components 
at high frequencies. In such a case we recommend to use a window function (e.g. the 
Bartlett window) while calculating the inverse transform: 



lL + J- m=-Lf3€C \ 



. e i(kf3-k)t„ 
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9c(k) = -f— £ £ 1 " T W {k *~ k)tm - (2-15) 

m=—L /3gC \ / 

If the distance between the eigenvalue and the limit of the interval on the real k— axis 

Ak 
is smaller than — , it gives rise to the spurious periodic behavior of the inverse 



Fourier transform. To avoid it one should calculate Fourier transformation with en- 
larged value of Ak, namely 



Ak -> Afc(l + a) 



(2.16) 



with a value for a ~ 0.05 . . . 0.15. 



As we have mentioned above the simplest and perhaps the most useful type of contour C 
is a rectangle. Let us define such a rectangle by its four corners (a — iS,b — i8, b + i5,a + i5) 
with < a < b and 5 > 0. Then the integrals Rel. (|2.12 ) can be somewhat simplified and 
thus the amount of numerical work reduced: 






\' a (x - 18) 

( x 

X a (x — iS) 



i5) n 



(2.17) 



Jo 2vr 



a=l 



\ a (b + iy) \ a {a + iy) 



As strange as it might sound, any of the known in the literature methods should encounter 
most troubles for problems with high symmetries, especially in higher than two dimensions, 
when high degeneracies are present. In such cases one has to determine not only the presence 
of an eigenstate, but also its degeneracy. The problems without symmetries are easier in 
this respect, as degeneracies are rare and accidental. However, once in a while an accidental 
degeneracy might appear. When tunneling occurs, very closely spaced doublets appear in 
the spectrum [12| . For the method presented here that does not seem to pose any challenges 
(except perhaps the actual resolution of such a doublet, but not the determination of its 
existence). The dimensionality of the problem appears to be largely irrelevant as well, except 
for the obvious increase of the dimensionality of the BOM. Some of these aspects will be 
illustrated in the next section. 

In many instances one is interested in evaluating the ground state energy of a many 
fermion system for a given shape and a given number of particles N, that is 



E. 
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2m 



S 2 (C), 



where the contour C has to be chosen so as to satisfy the condition 

N = JV(C). 



(2.18) 



(2.19) 



The contour could be a union of several contours, chosen in accordance with the "rules" 
discussed above. In order to evaluate the ground state energy E gs one needs the position of 
the Fermi level, which is determined by solving Eq. ( p,19| ). The exact position of the other 
eigenvalues are not needed and thus one has to determine exactly only a few eigenvalues 
around the Fermi level. 



III. A FEW PARTICULAR CASES 

The examples we are going to discuss here are particularly instructive. In many of 
these cases some, most or all of the calculations can be performed analytically. Thus these 
examples serve as an extremely useful guide through the potential numerical problems, some 
of which we have alluded to in the previous section. In particular, from the exact expressions 
for the BOM eigenvalues given below, it follows that for k — > ±ioo the integrand behaves 
as X' a (k)/X a (k) oc 1/k and thus it is not profitable to deform the contour C into a straight 
parallel to the imaginary £;-axis, as one might have been tempted to. 

A. 1 dimensional segment 

One can easily show that for a 1-dimensional segment of length L one can solve the 
problem exactly using the outlined method. The BOM has exactly two eigenvalues 

X 1 (k) = 2sin 2 ^, (3.1) 

A 2 (A;) = 2cos 2 — . 

The contour integrals ( |2.9[ ) and ( |2.12|) can be evaluated analytically and they lead to the 
expected exact results. One can see here explicitly that the roots of the BOM eigenvalues 
have double multiplicities, as we have stated in the previous section. 



B. Circular Billiard 

If one considers a circular billiard of radius R, the BOM elements can be given analytically 

0(6, 6') = J \2kR sin B -^- j . (3.2) 

Here Jo(x) is the cylindrical Bessel function of first kind and 9 and 6' determine the direction 
of two unit vectors in cylindrical coordinates in an obvious manner, k = (cos#, SU16 1 ). In a 
rather limited so far numerical implementation of the present method we have established 
that all of the exact solutions of the circular billiard are reproduced. 

One can show that for the circular billiard the exact BOM eigenvectors and BOM eigen- 
values (in the continuum limit when the number of basis plane waves is infinite and N — > oo) 
are c{9) = exp(im9)/\/2~K and \ m = J^(kR) respectively, with an arbitrary integer m. 
Therefore, the solution one obtains using the method described in the present paper leads 
to the exact result for the eigenvectors and eigenvalues of the Helmholtz equation with 
Dirichlet boundary conditions. Note again the fact that the roots of the BOM eigenvalues 
X m have double multiplicities. 



C. Annular Billiard 

We shall consider now the billiard of radius R with a circular hard core of radius a < R, 
with an offset 5 < R — a along the x-axis (see Fig. |l]). For the numerical implementation 
of the method for the annular billiard it is more convenient to use the basis wave function 



in the form [13 



i/> n (kr, 4> a ) = [Jm{kr)Y m (ka) - Y m (ka)J m (kr)} exp(imd a ), (3.3) 

where J m and Y m are cylindrical Bessel functions of the first and second kind, respectively. 
This implies that the boundary conditions on the inner circle are automatically fulfilled. 
One can convince oneself that the popular plane wave basis is inadequate in this case. 
The simplest way to see this is by attempting to use a plane wave basis to represent the 
present basis wave functions. If the Bessel functions of the first kind can be represented as a 
superposition of plane wave of a given energy, the Bessel function of the second kind cannot. 
Since the system is invariant with respect to the Ci symmetry the BOM matrix decomposes 
into two blocks with the following matrix elements: 

Omm> = 7T dd cos(m6 a )[J m (kr)Y m (ka) - Y m (ka)J m (kr)} x 

Z7T JO 

cos(m'<f) a )[J m i(kr)Y m >(ka) - Y m '(ka)J m >(kr)] 

^W = 7T d(j)S\n{m(j)a)[J m {kr)Y m {ka) - Y m (ka)J m (kr)} x 
Zn Jo 

sm(m'4> a )[J m/ (kr)Y ml (ka) -Y m >(ka)J m >(kr)], (3.4) 

cos 6 — 8 



where r = \J1 + 5 2 — 25 cos 0, cos 6 a = and the radius of the billiard is set equal 

r 
to 1 (see Fig. 0). 

The number of basis states needed for the calculations up to the some value k max 
can be estimated from the condition that the only states contributing to the solution 
should have at least one node associated with the radial motion inside the domain T> i.e. 
^r£{a,i+5)\4>n<N(k max ,r)\ 2 = 0. An increase of the basis affects the computation time. The 
main contribution comes from the diagonalization of the BOM matrix along the integration 
contour. In the table [I] we have shown the dependence of the computation time on the 
number of basis states ip n . Although the computation time is sensitive also to the precision 
of the integration both in the (r, 6)— and k— spaces, it depends approximately linearly on 
the number of integration points whereas it exhibits the quadratic behavior as a function of 
the number of basis states. 

In the Figure |2] we have shown results of calculations performed for a = 0.5 and 5 = 
0, 0.25. Since for the 2-dimensional system the density of states is in the first approximation 
independent of k 2 we have used the variable e = k 2 instead of k. The level density p was 
obtained using L = 8 Fourier components in each interval Ae = 60, see Eq. (|2.14|) . The 
number of particles and the energy are expressed in the form: 

N= rdep(e)= rdeJ2PcM 
Jo Jo t 

deg(e) = / de^# Ci (e), (3.5) 

Jo ~r 



E 



i-' 
o 



where /x is the chemical potential. Since we are rather interested in the fluctuating part of 
these quantities we have substracted the smooth behavior associated with the density po 
calculated from the Weyl formula: 

Po(e) = — A 1 ^. (3.6) 

Then the corresponding particle and energy fluctuations are given by: 

r »r/ \ nr/ x 1 — a 2 1 + a 
<JiV(e) = iV(e) - — ^e - -^ 

5E(N) = E(N)- r°dep (e)e, (3.7) 

Jo 

(3.8) 
where /i is determined by the condition: 



/-'(> 



iV= / dep (e). (3.9) 



o 



In Ref. |14j| . it was shown that the annular billiard is chaotic for S > 0. However the 



degree of chaoticity i.e. the fraction of the phase space occupied by chaotic trajectories 
depends on the eccentricity parameter 5. The system becomes fully chaotic for S = 1 — a 
when all the trajectories must hit the inner circle. Thus the shell effects visible in the Fig. 



for S = 0.25 are remnants of the ordered motion still existing in this case, see Refs. [15 



Associated orbits are characterized by the impact parameter L = sin a > a + 5 (see Ref. 



14| ) and forever encircle the inner disk. One can observe that the period of oscillations of 
5E is larger for 5 = 0.25. It is caused by the fact that the shortest periodic orbits giving 
rise to the shell effects at 5 = i.e. triangular and rectangular orbits, are destroyed in the 
5 = 0.25 case. 

Within the method presented in this paper we are also able to plot the density distribu- 
tion of particular states. In Fig. |3] we have presented the densities of two states for 5 = 0.0 
and S = 0.25 belonging to two representations of Ci denoted by (+) and (— ). 

D. n dimensional billiards 

For billiards in any dimensions the BOM elements have simple forms if either the en- 
tire boundary or only parts of it have spherical/ellipsoidal (or 2d-circular/ellipse) shape, 
planar/linear and/or any combinations of these shapes. The sphere/circle could be easily 
deformed into an ellipsoid as well and one would still get simple expressions for the BOM 
elements. In particular for a spherical 3d-cavity one can show, similarly to the case of 
the 2<i-circular cavity, that the BOM eigenstates are as expected the spherical harmonics, 
i.e. ci m (9,(j)) = Yi m (9,(j)), and the BOM eigenvalues are Xi m = jf(kR) for m = —I,..., I, 
with spherical Bessel functions instead of cylindrical ones. This again shows that the BOM 
method leads to the correct eigenstates. 



10 



IV. CONCLUSIONS 

The method introduced here appears flexible enough to allow essentially a foolproof cal- 
culation of the entire spectrum of an arbitrarily shaped billiard in any dimensions. We 
believe also that the method can be optimized as well so as to speed up numerical calcu- 
lations. In particular, the introduction of evanescent plane waves |L6j is one aspect which 
is worth a serious consideration, as it apparently leads to significant increase of the numer- 



ical accuracy [I0|. Other boundary conditions, besides Dirichlet, can be implemented in a 
straightforward manner too. We were in way rather surprised to realize that the present 
method (which of course is nothing but a consequence of Cauchy's theorem) is not appar- 
ently being used in numerical analysis, as its implementation is rather straightforward. The 
fact that this method is a "foolproof" one should make it a standard one whenever finding 
roots of a rather large class of one variable functions is somewhat of a challenge. 
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TABLES 

TABLE I. Computation time estimates for the annular billiard with 8 = 0.25. The contour C 
has been choosen to enclose the interval (960, 1030) on the real e axis. 

No. of basis states Time (sec.) -^(C) 

2 11.13 0.000 

4 15.38 0.000 

6 21.87 0.000 

8 49.49 0.000 

10 58.08 0.000 

12 67.44 1.001 

14 79.39 1.000 

16 90.82 2.001 

18 106.78 2.000 

20 121.18 1.999 

22 140.28 4.015 

24 159.56 5.000 

26 185.45 7.000 

28 206.25 7.001 

30 235.27 9.000 

32 258.33 9.001 

34 300.64 11.000 

36 319.14 11.000 

38 361.64 11.000 

40 393.90 11.000 



13 



FIGURES 
FIG. 1. Annular billiard. 

FIG. 2. The density of states (p), particle number fluctuations (N — N smoot h) and energy 
fluctuations (E — E smoo th) as a function of e = k 2 for annular billiard. 

FIG. 3. Selected densities of quantum states of the annular billiard (the sign in paranthenses 
refers to the C2 symmetry of the corresponding wave function): a) k = 26.1,5 = 0, (+); b) 
k = 26.1,(5 = 0, (-); c) k = 25.354,5 = 0.25, (+); d) k = 25.096,5 = 0.25, (-). 
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